{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal,ogr,osr\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算水体的栅格像元数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pixel_Number(tif_path):\n",
    "    dataset = gdal.Open(tif_path)  #输入提取完成水体后的二值栅格\n",
    "    band = dataset.GetRasterBand(1)\n",
    "    array = band.ReadAsArray()\n",
    "    lake_area = np.sum(array == 1)\n",
    "    print(\"湖泊面积(像元数)：\", lake_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 将水体区域二值栅格转为矢量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Raster_to_Polygon(RasterPath, VectorPath):\n",
    "    ds_raster = gdal.Open(RasterPath)                        #读取栅格数据\n",
    "    band_raster = ds_raster.GetRasterBand(1)                 #获取需要转为矢量的波段\n",
    "    proj_raster = osr.SpatialReference()                     #获取栅格数据的投影信息\n",
    "    proj_raster.ImportFromWkt(ds_raster.GetProjection())     #将栅格数据的投影信息赋值给矢量\n",
    "    driver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "    if os.path.exists(VectorPath):                           #若文件已经存在，就删了重做一遍\n",
    "        driver.DeleteDataSource(VectorPath)\n",
    "    polygon = driver.CreateDataSource(VectorPath)            #创建数据资源\n",
    "    layer_polygon = polygon.CreateLayer(\"Shp\", srs=proj_raster, geom_type = ogr.wkbMultiPolygon)  #创建图层，定义多面\n",
    "    new_field = ogr.FieldDefn('value', ogr.OFTReal)          #给目标shp文件添加一个字段，用来存储原始栅格的像元值\n",
    "    layer_polygon.CreateField(new_field)\n",
    "    gdal.FPolygonize(band_raster, None, layer_polygon, 0)    #执行栅格转矢量操作\n",
    "    polygon.SyncToDisk()\n",
    "    polygon = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 计算面积"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Get_Polygon_Area(path_shp):                    #传入水体边界的矢量数据\n",
    "    driver = ogr.GetDriverByName(\"ESRI Shapefile\") #创建数据驱动\n",
    "    ds = driver.Open(path_shp, 1)                  #创建数据资源\n",
    "    layer = ds.GetLayer()\n",
    "    new_field = ogr.FieldDefn(\"Area\", ogr.OFTReal) #创建新的字段\n",
    "    new_field.SetWidth(32)                         #设置字段宽度（存储长度）\n",
    "    new_field.SetPrecision(16)                     #设置字段精度\n",
    "    layer.CreateField(new_field)\n",
    "    number = 0\n",
    "    Total_Area = 0\n",
    "    for feature in layer:\n",
    "        geom = feature.GetGeometryRef()\n",
    "        geom2 = geom.Clone()\n",
    "        area = geom2.GetArea()                     #默认为平方米\n",
    "        area = area / 1000000                      #将单位转为平方千米\n",
    "        feature.SetField(\"Area\", area)\n",
    "        Total_Area += area                         #计算水体总面积\n",
    "        layer.SetFeature(feature)\n",
    "    ds = None\n",
    "    print(f\"水体总面积：{Total_Area}平方千米\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 调用"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'gdal' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32md:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\Code\\FinalTest_tmp\\Area (1)\\area.ipynb 单元格 9\u001b[0m line \u001b[0;36m3\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m TIF_PATH\u001b[39m=\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mD:\u001b[39m\u001b[39m\\\u001b[39m\u001b[39mLenovo\u001b[39m\u001b[39m\\\u001b[39m\u001b[39mDesktop\u001b[39m\u001b[39m\\\u001b[39m\u001b[39m云南大学\u001b[39m\u001b[39m\\\u001b[39m\u001b[39m开源地理信息系统\u001b[39m\u001b[39m\\\u001b[39m\u001b[39m期末实验数据\u001b[39m\u001b[39m\\\\\u001b[39;00m\u001b[39mtmp\u001b[39m\u001b[39m\\\\\u001b[39;00m\u001b[39mtest7.TIF\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m Vector_Path\u001b[39m=\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mD:\u001b[39m\u001b[39m\\\u001b[39m\u001b[39mLenovo\u001b[39m\u001b[39m\\\u001b[39m\u001b[39mDesktop\u001b[39m\u001b[39m\\\u001b[39m\u001b[39m云南大学\u001b[39m\u001b[39m\\\u001b[39m\u001b[39m开源地理信息系统\u001b[39m\u001b[39m\\\u001b[39m\u001b[39m期末实验数据\u001b[39m\u001b[39m\\\\\u001b[39;00m\u001b[39mtmp\u001b[39m\u001b[39m\\\\\u001b[39;00m\u001b[39mwater1.shp\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m----> <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m Pixel_Number(TIF_PATH)\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m Raster_to_Polygon(TIF_PATH,Vector_Path)\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=4'>5</a>\u001b[0m Get_Polygon_Area(Vector_Path)\n",
      "\u001b[1;32md:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\Code\\FinalTest_tmp\\Area (1)\\area.ipynb 单元格 9\u001b[0m line \u001b[0;36m2\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mPixel_Number\u001b[39m(tif_path):\n\u001b[1;32m----> <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m     dataset \u001b[39m=\u001b[39m gdal\u001b[39m.\u001b[39mOpen(tif_path)  \u001b[39m#输入提取完成水体后的二值栅格\u001b[39;00m\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m     band \u001b[39m=\u001b[39m dataset\u001b[39m.\u001b[39mGetRasterBand(\u001b[39m1\u001b[39m)\n\u001b[0;32m      <a href='vscode-notebook-cell:/d%3A/Lenovo/Desktop/%E4%BA%91%E5%8D%97%E5%A4%A7%E5%AD%A6/%E5%BC%80%E6%BA%90%E5%9C%B0%E7%90%86%E4%BF%A1%E6%81%AF%E7%B3%BB%E7%BB%9F/Code/FinalTest_tmp/Area%20%281%29/area.ipynb#X11sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m     array \u001b[39m=\u001b[39m band\u001b[39m.\u001b[39mReadAsArray()\n",
      "\u001b[1;31mNameError\u001b[0m: name 'gdal' is not defined"
     ]
    }
   ],
   "source": [
    "TIF_PATH=\"D:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\期末实验数据\\\\tmp\\\\test7.TIF\"\n",
    "Vector_Path=\"D:\\Lenovo\\Desktop\\云南大学\\开源地理信息系统\\期末实验数据\\\\tmp\\\\water1.shp\"\n",
    "Pixel_Number(TIF_PATH)\n",
    "Raster_to_Polygon(TIF_PATH,Vector_Path)\n",
    "Get_Polygon_Area(Vector_Path)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv-2",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
